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We present a perturbative treatment of the response properties of insulating crystals under a dc 
bias field, and use this to study the effects of such bias fields on the Born effective charge tensor and 
dielectric tensor of insulators. We start out by expanding a variational field-dependent total-energy 
functional with respect to the electric field within the framework of density-functional perturbation 
theory. The second-order term in the expansion of the total energy is then minimized with respect 
to the first-order wave functions, from which the Born effective charge tensor and dielectric tensor 
are easily computed. We demonstrate an implementation of the method and perform illustrative 
calculations for the III-V semiconductors AlAs and GaAs under finite bias field. 

PACS numbers: 71.15.-m, 71.55.Eq, 77.22.-d, 78.20.Ci. 
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I. INTRODUCTION 

The dielectric tensor and Born (or dynamical) effec- 
tive charge tensor are of fundamental importance in un- 
derstanding and modeling the response of an insulator 
to an electric field. 1 They give, respectively, the first- 
order polarization and atomic force appearing in response 
to a first-order change in the macroscopic electric field. 
While one is most often interested in evaluating these 
response tensors at zero field, there is increasing inter- 
est in finite-field properties. For example, the study 
of bulk ferroelectrics 2-4 and of ferroelectric films 5 and 
superlattices 4 ' 6 in finite field, and of lattice vibrations in 
polar crystals in finite field, 7 have recently generated in- 
terest. While it may sometimes be reasonable to model 
the dielectric behavior by assuming that the dielectric 
and Born effective charge tensors have a negligible depen- 
dence on the bias field, it is important to be able to quan- 
tify such approximations and to compute the field depen- 
dence when it is physically important to do so (e.g., for 
describing non-linear optical phenomena such as second- 
harmonic generation) . 

Density-functional perturbation theory (DFPT) 8,9 
provides a powerful tool for calculating the second-order 
derivatives of the total energy of a periodic solid with re- 
spect to external perturbations such as atomic sublattice 
displacements or a homogeneous electric field. In con- 
trast to the case of sublattice displacements, for which 
the perturbing potential remains periodic, the treatment 
of homogeneous electric fields is subtle because the cor- 
responding potential acquires a term that is linear in real 
space, thereby breaking the translational symmetry and 
violating the conditions of Bloch's theorem. For this rea- 
son, electric-field perturbations have often been studied 
in the past using the long-wave method, in which the 
linear potential resulting from the applied electric field 
is obtained by considering a sinusoidal potential in the 
limit that its wave vector goes to zero. In this approach, 
however, the response tensor can only evaluated at zero 
electric field, and it also requires as an ingredient the 
calculation of the derivatives of the ground-state wave 



functions with respect to wave vector. 

Recently, Nunes and Gonze introduced an electric- 
field-dependent energy functional expressed in terms of 
the Berry-phase polarization. 10 This approach was ini- 
tially introduced in order to provide an alternative frame- 
work for the DFPT treatment of electric-field perturba- 
tions (evaluated at zero field) in which the long-wave 
method is entirely avoided. More recently, it has been 
pointed out that the Nunes-Gonze functional could also 
serve as the basis for a calculation of the "ground- 
state" properties in finite electric field. 11,12 (Here, the 
phrase "ground state" is used advisedly; because of Zener 
tunneling, the state of interest is actually a long-lived 
resonance. 13 ) In this approach, the energy functional is 
minimized with respect to a set of field-polarized Bloch 
functions that form a natural representation of the one- 
particle density matrix even though they are no longer 
eigenstates of the Hamiltonian. n ' 13 The introduction of 
this approach has also made possible the calculation of 
the dielectric and Born effective charge tensors at finite 
electric fields using finite-difference methods. 11 

In a recent paper 7 we developed a perturbative 
method, within this framework, for computing the 
phonon properties of insulators at finite electric fields. 
The starting point was the Nunes-Gonze electric-field- 
dependent energy functional, which represents the effect 
of the electric field by including its coupling to the Berry- 
phase polarization. 10 This total-energy functional was ex- 
panded up to second order in atomic displacements. The 
linear response of the field-polarized Bloch functions to 
the atomic displacements was obtained by minimizing the 
second-order term in the expansion of the total-energy 
functional with respect to the first-order changes in the 
Bloch functions. Finally, the force-constant matrix was 
constructed based on these first-order Bloch functions. 
This method provides a tractable and efficient computa- 
tional scheme for computing phonon properties at finite 
electric field, and suggests that a similar treatment of 
other response properties of insulators in finite electric 
field should be possible. 

In this paper, we follow the approach of Ref. 7 to de- 
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velop a method for computing the dielectric and Born 
effective charge tensors at finite electric field. Again us- 
ing the Nunes-Gonze energy functional, 10 we compute 
the first-order responses of the electronic wave functions 
to a small change in the electric field. We then use these 
to construct the second-order derivatives of the total en- 
ergy with respect to electric field (to give the dielectric 
tensors evaluated at non-zero field) and the mixed deriva- 
tives with respect to electric field and atomic sublatticc 
displacement (to give the Born effective charge tensors 
evaluated at non-zero field). 

The paper is organized as follows. In Sec. II, the 
second-order perturbation expansion of the total energy 
functional with respect to electric fields is derived and 
the steepest-descent directions are identified. The ex- 
pressions for computing the dielectric and Born effec- 
tive charge tensors are also given. In Sec. Ill, we de- 
scribe the implementation of the approach in the ABINIT 
code package, 14 and present test calculations for the III- 
V semiconductors AlAs and GaAs. (Since we are mainly 
interested in the purely electronic effects here, we do not 
include the strains or sublattice displacements that might 
occur in response to the electric field; these could easily 
be included by employing structural-relaxation methods 
at finite field. 11 ) By comparing with the results of finite- 
differences calculations, we demonstrate the correctness 
of the new formulation and the internal consistency of the 
theory. A brief summary and conclusion are presented in 
Sec. IV. 



II. METHOD 

A. Perturbation expansion of the enthalpy 
functional 

We start from the electric enthalpy functional 10,11 

F[R;V;£]=£ K s[R;V>]-^-P[V>] . W 

where R, £, £1 and P are, respectively, the atomic coor- 
dinates, the electric field, the cell volume, and the macro- 
scopic polarization, Fks is Kohn-Sham energy functional 
at zero electric field, and atomic units are used through- 
out. After minimizing this functional, the field-polarized 
Bloch functions tp may be regarded as depending implic- 
itly on the electric field £. Our treatment of this func- 
tional will parallel the treatment given in our previous 
Ref. 7. 

In the present case, we take the electric field £ to con- 
sist of two parts, a finite part £(°> and a small variation 
S£. In the following, we consider the perturbation ex- 
pansion of the functional in Eq. (1) with respect to the 
small variation S£ under the orthonormality constraints 

{lpmk\lpnk) = $mn ■ (2) 

The wave functions are to be relaxed, subject to these 
constraints, in such a way as to minimize the electric 



enthalpy functional 

F = F KS + F BP + f LM , (3) 

where Fks = Fks is the Kohn-Sham energy (as it would 
be calculated at £ = 0), Fbp = —Of • P contains the 
coupling of the Berry-phase polarization P to the electric 
field, and the constraint is implemented by the inclusion 
of the Lagrange-multiplier term J^lm ■ The first and last 
of these terms are given by 

n occ 

fextlVVik) + F Hxc [n] (4) 

^ kn 

and 

j occ 

FhM = -jy- ^ A k ,mn((?/'<nk|V'nk) - S mn ) (5) 
k,mn 

where / is the spin degeneracy (normally /=2), Nk is the 
number of fc-points, and Ak, m « is the matrix of Lagrange 
multipliers. In a notation similar to that of Ref. 7, the 
second term may be written as 

i=l JV _L k 

Here are the three primitive real-space lattice vectors, 
and the mesh of Nk fc-points is defined by mesh vectors 
gi = bj/iVW where is the reciprocal lattice vector 
dual to a;. Thus, N k = N^N^N^, and we also define 
= Nk/N^ as the number of fc-point strings running 
in direction i. Finally, 

-CW = Im In det SW (?) 
where the overlap matrix is defined as 

(SW)nm = (Umk|«nk'} • (8) 

In order to obtain the desired response properties, we 
now wish to expand the finite-field enthalpy functional 
Fks up to second order in the electric field. We shall as- 
sume for the moment that the electric field is applied in 
Cartesian direction a only. The expansion of Fks with 
respect to atomic displacements was already obtained in 
Ref. 7, and the expansion with respect to electric field 
can be carried through in a very similar way. Indeed, the 
second-order expansions of Fks and Flm can essentially 
be transcribed from Ref. 7 with the first-order wave func- 
tions with respect to displacement replaced here by the 
first-order wave functions with respect to electric field, 
giving 

p( 2 ) _ 1 d 2 F KS 
KS " 2 d£l 

j, OCC 

= w EE<*4i T+t; exti<k> 

k k n=l 

(9) 
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and 
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p(2) 
LM 
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( W nkl M nk) ■ 



(10) 



As in Ref. 7, terms that can be eliminated by use of the 
"2n+l theorem" (e.g., (u^ £a \T + v ext \u^l)) have been 
dropped. The the first-order wave functions are 



) = 



d£„ 



and the second-order En KC are 

9 2 £'hxc 



E 



£ a £ a 
Hxc 



2d£ a d£ a 



(11) 



(12) 



In these and subsequent equations, the partial deriva- 
tives indicate that the structural coordinates R are being 
held fixed (while, however, the wave functions |u„ k ) arc 
allowed to vary). 

The second-order expansion of Fbp with respect to 
electric field requires somewhat more care. We find 



7 (2) 

BP 



1 d 2 F BP 

2 d£l 

n d 2 (£-p) 



d£l 



= -0(e o 



+ £(0). P £«£«) } (13) 



where e a is the unit vector along Cartesian direction a. 
The first term in the last line of Eq. (13) is special to the 
case of the electric-field perturbation, while the second 
term can be derived in close correspondence to the case 
of displacement perturbations in Ref. 7. The first-order 
variation of P with field £ a is 



- e/ \ - a, (i) 
2ttO ^ ^ 



_a 

k 

and its second-order variation is 
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k,k+ gl 



y£ a £ a 



e J_ sr V r>( 2 ) 

i=l JV ± k 



where 



^k^+g.Qk+g^k 



(14) 



(15) 



(16) 



and 



(2) 

k,k+ gi 



= ImTr 



2 ^k,k+ g , ( 3k+g 1 ,k 



5 'k,k+g I ( 3k+g 1 ,kS' k!k+ g.Qk+g i ,k • (17) 



In these equations, 'Tr' indicates a trace of the bracketed 
matrix over band indices, and Q, and are de- 

fined with respect to the series expansion of the overlap 
matrix via 



c2 c(2) 
C Q J kk' 



(18) 



Qkk< = [sgU- 1 



(19) 



The first- and second-order expansions of the overlap ma- 
trix take the form 



4,k',mn — ( U mk\ U nl') + ( U lnL\ U nk') ( 2 0) 



and 



'mkrnk'/ • (21) 

In the last equation above, terms like (uf"£ a |«^° k /) have 
again been dropped by virtue of the "2n + 1 theorem." 



First-order wave functions with respect to 
electric-field perturbation 

The second-order term in the expansion of the energy 

F (2) of 

LM oi 



(2) 



BP 



functional, given by the sum — F^g 
the expressions in Eqs. (9), (13), and (10) respectively, is 
minimized with respect to the first-order wave functions 
\ u nk) usm g standard conjugate-gradient methods. The 
steepest-descent direction is obtained from the gradient 



of F^ with respect to 



the form 



Sikh 



whose contributions take 



SF 



(2) 



KS 



Su 



nk 



J_ 



Su 



. (22) 



5F™ ie/A ^ 



nk 



-i=1 iv J_ 
. , 3 „ 

+ 4^ Pmk.k+gi) - l^mk.k-gj J 



and 



Here 



^LM _ / (°)| 7 A\ 



(23) 



(24) 



CW = ( |« k ?)Q k , k - kk'^k'k^kk'Ok'k ) ■ (25) 



Z>mkk' = (lu^^Ok'k) » 

V / m 



(26) 



and is the diagonal zero-order matrix of Lagrange 
multipliers. Convergence of the conjugate-gradient pro- 
cedure yields a set of first-order wave functions |u nk ). 
These then become the essential ingredients for con- 
structing the dielectric and Born charge tensors as dis- 
cussed below. 
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C. Dielectric permittivity tensor 



The dielectric permittivity tensor can be written as 



Zap = <W + 47TXa/3 



(27) 



where the electric susceptibility tensor Xa/3 at a finite 
electric field is defined as 



XaO 



1 d 2 F(£ ) 



£1 d£ a d£fj 
dP c 



£=£(0) 



= e G 



£=£(0) 



(28) 



The derivative P £f3 of the polarization with respect to 
electric field is already given by Eq. (14). Since the first- 
order wave functions \uf^) have already been obtained in 
Sec. II B, it is straightforward to evaluate Eq. (28) and 
thus obtain the polarizability and permittivity. 

The dielectric responses above are the static responses 
computed with atomic coordinates frozen. That is, they 
correspond to the dielectric response that would be mea- 
sured at frequencies low compared to electronic frequen- 
cies but high compared to any infrared-active phonon 
modes. The true static susceptibility could be computed 
by including the lattice displacements (and, if appropri- 
ate, the piezoelectric strains) using, e.g., the methods of 
Ref. 15. 



Alternatively, Eq. (29) can be computed as the deriva- 
tive of the polarization with respect to the displacement, 



dP a 



(33) 



Here P Tk ^ takes a form very similar to that of Eq. (14), 



except that the first-order changes 



in the wave 



functions in response to an electric field are replaced by 



the corresponding changes 



in response to a sub- 
3 ) has 



lattice displacement. The computation of the \u 
already been described in detail in Rcf. 7. 

The computation of the first-order derivatives of the 
wave functions is typically the most time-consuming step 
of the linear-response calculation. Therefore, for a com- 
plicated unit cell with many atoms M per cell, the com- 
putation of the three derivatives \u £a ) will be much 
cheaper than that of the 3M derivatives Im 1 ""'* 9 ), and the 
method of Eq. (32) will therefore be significantly faster 
than the method of Eq. (33). In the special case that 
the displacement derivatives \u T "- 13 } have already been 
computed for some other reason (e.g., for the purpose of 
computing the phonon frequencies in finite field), the use 
of the latter method may be advantageous. In any case, 
a comparison of the two methods should provide a useful 
check on the internal consistency of the theory and its 
computational implementation. 



D. Born effective charge tensor 

The electronic contribution to the Born effective charge 
tensor at finite electric field takes the form 



d 2 F(£) 
^ ~ d£ a dr K . p 



£=£(«) 



(29) 



This expression can be calculated equivalently in two 
different ways. First, introducing the force f KM = 
— dF{£) / dr K>a acting on atom k in direction a, it can 
be written as 



Z. 



K,a{3 



d£ a 



(30) 



Using the Hellmann-Feynman theorem, the expression 
for the force is given as 

» occ 



k n=l 



and taking an additional derivative with respect to elec- 
tric field yields 

_ , occ 

z t«? = ^EE< u iki( T +^r^i< k > • (32) 



k ra=l 



This has essentially the same form as Eq. (43) in Ref. 9, 
except that here the zero-order wave functions are al- 
ready polarized by the preexisting finite electric field. 



III. TEST CALCULATIONS FOR III-V 
SEMICONDUCTORS 

In order to check our method, we have performed 
test calculations on two prototypical III-V semiconduc- 
tors, AlAs and GaAs, for which the electronic contri- 
bution to the polarization is typically comparable to 
the ionic contribution. 4 The calculation is carried out 
using the planewave-pseudopotential method based on 
density-functional theory with local-density approxima- 
tion (LDA). We use Troullier-Martins norm-conserving 
pseudopotentials 16 in which the 3d states on the Ga and 
As atoms are treated as core states. (The omission of 
the semicore 3d states from the valence on the Ga atom 
may limit the accuracy of the Ga pseudopotcntial some- 
what.) A 16 x 16 x 16 Monkhorst-Pack mesh is used for 
the fc-point sampling. More computational details can be 
found in our preceding paper. 7 

The calculation of the dielectric permittivity tensor 
and the Born effective charge tensor is carried out in three 
steps. First, a ground-state calculation at finite electric 
field is performed using the Berry-phase approach 11 im- 
plemented in the ABINIT code, and the field-polarized 
Bloch functions are stored for the later linear response 
calculation. Second, the linear response calculation is 
carried out to obtain the first-order response of Bloch 
functions. Third, the matrix elements of the dielectric 
and Born effective charge tensors are computed using 
these first-order responses. 
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TABLE I: Calculated electronic dielectric constants of AlAs 
and GaAs at zero field, and changes resulting from an electric 
field of 3.08 x fO 8 V/m along the [fOO] direction. 'LR' and 
'FD' denote the results of linear-response [Eq. (28)] and finite- 
difference calculations, respectively. 





too 


Aeoo,23 


Aeoo,n 




AlAs LR 


9.681 


0.039 


0.027 


0.013 


FD 


9.681 


0.040 


0.027 


0.013 


GaAs LR 


13.315 


0.202 


0.211 


0.104 


FD 


13.319 


0.203 


0.207 


0.098 



TABLE II: Calculated cation Born effective charges of AlAs 
and GaAs at zero field, and changes resulting from an electric 
field of 3.08 x 10 s V/m along the [100] direction. 'LR' and 
'FD' denote the results of linear-response [Eq. (32)] and finite- 
difference calculations, respectively. 







Z* 


AZ 2 * 3 
(xlO -3 ) 


AZ* n 
(xl0~ 3 ) 


AZ 3 *3 

(xlO -3 ) 


AlAs 


LR 


2.110 


17.23 


-0.06 


-0.13 




FD 


2.110 


17.22 


-0.05 


-0.11 


GaAs 


LR 


2.186 


52.88 


-3.42 


-3.17 




FD 


2.186 


52.83 


-3.36 


-3.14 



The first column of Table I shows the calculated elec- 
tronic dielectric constants of AlAs and GaAs at zero 
electric field, and the remaining ones show the nonzero 
changes in the dielectric tensor elements after the appli- 
cation of an electric field £(°) of 3.08 x 10 8 V/m along 
the [100] direction. The results obtained with the linear- 
response approach of Eq. (28) are compared with those 
calculated by finite differences. In the latter case, polar- 
izations are computed at several values of the electric field 
in steps of 3.08 x 10 5 V/m, and the dielectric tensor is cal- 
culated using a finite-difference version of Eq. (28). It can 
be seen that the agreement between the linear-response 
and the finite-difference results is excellent, demonstrat- 
ing the internal consistency between the two approaches. 

In Table II we present similar results for the cation 
Born effective charges of the same two materials, first at 
zero field and then again under application of a field of 
of 3.08 x 10 8 V/m along the [100] direction. The 
linear-response results were obtained using Eq. (32), but 
we also computed the corresponding values using Eq. (33) 
and found agreement between the two linear-response ap- 
proaches with a maximum fractional error smaller than 
10 -6 for all values reported. For the finite-difference com- 
parison, the polarizations were computed at several val- 
ues of the atomic displacements in steps of 10~ 3 Bohr and 
the Born charge tensors were calculated using a finite- 
difference version of Eq. (33). It can again be seen the 
agreement between the linear-response and the finite- 
difference results is excellent. 

We emphasize that the values of Aeoo and AZ* re- 
ported in Tables I and II are purely electronic or "frozen- 
ion" ones - that is, the sublattice displacements that 
would be induced by a truly static electric field £(°) are 
not included. 

The values of eoo and Z* reported in Tables I and II 
arc in good agreement with other theoretical values in 
the literature 17-19 and with experiment. The symmetry 
is such that the applied electric field along x breaks the 
degeneracy between the diagonal elements of the £oc and 
Z* tensors so that 6oo.ii ^ £00,22 = £00,33 and Z* lx 7^ 
Z22 = -Z33, and introduces non-zero off-diagonal elements 
£00,23 = £oo, 32 and Z% 3 = Z| 2 . 

and 

00,11? 



Symmetry considerations also imply that e 



00,23 

Z23 should appear to first order in £(°\ while Ae 



Aeoo,33, AZ* l7 and AZg* 3 should be quadratic in £^°\ 
This is confirmed by our numerical calculations. Indeed, 
by repeating calculations like those shown in Tables I and 
II for several values of £^ and fitting to obtain the coeffi- 
cients of the linear and quadratic dependence, we can ex- 
tract information about the nonlinear dielectric response 
and the Raman tensor. The second- and third-order non- 
linear dielectric tensors are defined as 



1 



X123 — o 



2 88x883 



and 



(3) 

Xiiu - g 



1 8 3 P 1 
~8Ej 



1 <9X23 

2 8£i 



1 d 2 X ii 
6 



(34) 



d£{ 



(35) 



while the Raman polarizability tensor is defined by 

8Z 23 



d 2 h 



881883 8 81 



(36) 



where f is the force on the cation sublattice induced 
by the electric field. In practice, we calculate X23, 
X11, and Z23 for a scries of finite electric fields oriented 
along the 21-axis with values of £(°) ranging from zero 
to 5.14 x 10 8 V/m in increments of one-fifth of the max- 
imum value. Fitting these data to a polynomial in £(°) 
then gives the values of X1231 Xnin and cuto- Note that 
axo can alternatively be expressed as 



cvto 



n 



dx 



23 



8t\ 



(37) 



where T\ is a cation sublattice displacement and X23 is 
computed at zero field. We have also computed cvto 
by fitting to a series of calculations of this type, and 
find values of ctxo that agree with those obtained from 
Eq. (36) within 0.3%. 

(2) 

The results for the X123 and cuto values as computed 
from Eqs. (34) and (36) are presented in Table III for 
AlAs, together with some previous theoretical and exper- 
imental values for comparison. In view of the fact that 
the calculation of higher-order tensor elements tends to 
be delicate, the agreement is generally quite good. In 
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TABLE III: Values of second-order dielectric susceptibility 
and Raman matrix elements in AlAs, as defined by Eqs. (34) 
and (36) respectively, compared with previous theory and ex- 
periment. 



X [% (pm/V) 



|qto| (A 2 ) 



Present work 
Theory," Ref. 11 
Theory, 6 Ref. 20 
Theory," Ref. 21 
Theory, 6 Ref. 22 
Experiment, Ref 23 



62 
64 
70 
79 

78±20 



8.0 

8.5 
9.0 
7.4 



"Using finite-difference approach. 
6 Using (2n + l)-theorem approach. 



particular, Veithen et al. 20 have shown (see their Fig. 1) 

(2) 

that the results for X123 can t> e quite sensitive to the 
method of discretization in /c-space and the fineness of the 
fc-point mesh. For GaAs we find X123 = 293 pm/V and 
q:to = —24.1 A 2 (which is close to the value in Ref. 22), 
but these numbers are of questionable accuracy because 
of our use of a Ga pseudopotcntial that does not include 
the 3d semicore orbitals in the valence. We obtain X1111 
values of 3.90 and 33. 8x 1CT 11 esu for AlAs and GaAs, 
respectively. We are not aware of previous theoretical 
values of X1111 with which to compare; this quantity is 
beyond the reach of the "2n + 1" theorem using first- 
order wave function responses only, and so is difficult 
to compute by pure DFPT methods. Experimental val- 
ues ranging from 3.9 to 18xl0~ n esu for GaAs 24 can be 
found in the literature. 

The discrepancies noted above between theory and the- 
ory, and between theory and experiment, may have many 
possible causes. In addition to some of the computa- 
tional and convergence issues mentioned above, the ad- 
equacy of the LDA approximation itself is also a seri- 
ous question. Because the LDA tends to underestimate 
gaps, some authors have included a so-called "scissors 
correction" in order to widen the gap artificially; this 
tends to decrease the magnitude of response tensors. 25 
On the experimental side, the difficulty in obtaining re- 
producible results is surely also an issue. Nevertheless, 
we emphasize that the relative accuracy of the values re- 
ported in Tables I and II, which were done under the 



same computational conditions (same pseudopotentials, 
fc-point meshes, etc.), demonstrates the correctness of our 
new finite-field linear-response formulation and the inter- 
nal consistency of the computational framework that we 
employ. 



IV. SUMMARY 

We have developed a linear-response method for com- 
puting dielectric constants and Born effective charges in 
the presence of a finite electric field. We have demon- 
strated the reliability of our approach by implementing 
it in the context of the ABINIT code package 14 and per- 
forming test calculations on two III-V semiconductors, 
AlAs and GaAs. We have confirmed that the results cal- 
culated using the new linear-response approach are con- 
sistent with those obtained from finite-difference calcula- 
tions carried out within the same framework. In general, 
our results are also in good agreement with other theo- 
retical calculations and with experiment. 

A major advantage of the present approach is that, un- 
like the conventional long-wave linear-response method, 8 
it can be applied to obtain response tensors in finite elec- 
tric field. While it is possible to obtain similar infor- 
mation from a set of finite-difference calculations car- 
ried out for some chosen set of applied electric fields, the 
linear-response approach is more direct, and it avoids 
the troublesome truncation errors that may arise in a 
finite-difference approach. In the future, it may be of 
interest to extend the finite-field DFPT treatment not 
just to phonon perturbations (presented in Ref. 7) and 
electric-field perturbations (presented here), but also to 
other perturbations such as those associated with strain 
or chemical composition. Taken together, these develop- 
ments should allow for much greater flexibility in the cal- 
culation of materials properties of insulators under elec- 
trical bias and facilitate the study of higher-order non- 
linear dielectric properties. 
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